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ABSTRACT 

We present Smoothed Particle Hydrodynamics (SPH) simulations of molecular cloud 
formation in spiral galaxies. These simulations model the response of a non-self- 
gravitating gaseous disk to a galactic potential. The spiral shock induces high densities 
in the gas, and considerable structure in the spiral arms, which we identify as molec- 
ular clouds. We regard the formation of these structures as due to the dynamics of 
clumpy shocks, which perturb the flow of gas through the spiral arms. In addition, the 
spiral shocks induce a large velocity dispersion in the spiral arms, comparable with the 
magnitude of the velocity dispersion observed in molecular clouds. We estimate the 
formation of molecular hydrogen, by post-processing our results and assuming the gas 
is isothermal. Provided the gas is cold (T ^ 100 K), the gas is compressed sufficiently 
in the spiral shock for molecular hydrogen formation to occur in the dense spiral arm 
clumps. These molecular clouds are largely confined to the spiral arms, since most 
molecular gas is photodissociated to atomic hydrogen upon leaving the arms. 
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1 INTRODUCTION 

Molecula r clouds are the sites of star formation in spiral 
galaxies (lReddishl[l975l : IShu et all [l987l: iBlitz Williamsl 
° I1999I I. As such, they play an important role in the gov- 
erning the properties of galaxies and their evolution. The 
properties of molecular clouds have been extensively stud- 
ied in term s of their emission, k inematics and interna l struc- 
tur es (e.g. Dame et al. "198^; 'Solomo n et allll98/t iBruntI 
1200 3: Wilson et al. 2005). Unfortunately they are still poorly 
understood in terms of their formation and evolution. 
For example, it is unclear whether the molecular gas is 
formed in situ from atomic gas or whether molecular clouds 
are due to the coagulation of pre-existing molecular gas 
(iBlitz R.oso lowskv 2004; Prinde et al. 2001). Another im- 
portant question is whether molecular clo uds are long lived 
or just very inefficient at forming stars (|Elmegreenl bOOCt 
iHartmann et al] 1200 J) . Inefficient star formation can be 
potentially explained if the clouds are not globally self- 
gravitating (Clar k BonnelJ l2004l : IClark et all l200.tl . In 
this case, locally bound regions can explain the star for- 
mation while the vast majority of the gas escapes the cloud 
as it is unbound. 

Molecular clouds are found to be primarily located 
in spiral arms of galaxies with an arm to interarm ratio 
varying from a few in the inner regions of some galax- 
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ies (IVogel et al.ll98^lAdler WestDfahil996l :lR^ 
ll996l;"Brouillet et al."l998') to greater than 20 in th 



leuter et al.l 
L the oute r re- 

gk)ns of the Milky Way (Heyer k Terebey 1998; Digel e tHI 
119961) . Since spiral arms are generally defined by the pres- 
ence of star formation, this concentration of GMCs (gi- 
ant molecular clouds) to the spiral arms in the Galaxy 
is not surprising. The total fraction of gas in molecular 
form also varies from 2 per cent in M33 (Engargiola et alj 
[200j) to 75 per cent or more in M51 (fScoville Youna 
"l983f;|G arcia- Burillo et al.1ll993h . However, the detection of 
CO may be somewhat underestimated due to the uncer- 
tainty in th e ratio of CO(l- O) luminosity to H2 mass (the 
'X factor') ([Kaufman et al.l 1999). Interestingly, the mass 
spectra of GMCs in galaxies are ver y similar with dN oc 
m~^dm, 7 in the range of 1.5 to 1.8 (ISolomon et al.lll987l : 
iHever Terebevl [r998V Typical molecular cloud masses i n 
the Milky Way are 100-1000 M© (Hever Terebevlll998^ . 



have been obs erved in the Milky 
and M51 (iRand KulkarnJ 



thoug h GMCs up to 10^ M 
Way ([Solomon et al.l 

There are several possibiliti es regarding the f ormation 
of giant molecular clouds (see lElmegreenI Il990l and ref- 
erences therein). Where the content of molecular gas is 
high, e.g. M51 or the centres of galaxies, GMCs could 
form throu gh the agglomeration o f pre-existing molecular 
clouds (^e.g. IScoville Hershlll979l: ICowidll98d iTomisakal 
Il984l : iRoberts Stewart! 119871) . In this scenario, the ISM 
may be predominantly molecular, but GMCs are only 
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observed where the gas has coalesced into larger struc- 
tures heated by star formation. However this is un- 
likely to occur in the outskirts of galaxies or for exam- 
ple, M33, where l ittle molecular gas is thought to exist 
aside from GMCs (jEngargiola et alJl^ll^ . Then other pro- 
cesses, most commonly gravit ational /magneti c instabilities 
(Elmegreen 1979, 1982; Balbu s ^, Cowi3ll98,4 iBalbus 1988; 
|Elrneg^en||1989, 1991, 1994) or compression from shocks 
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j; iKoyama InutsukaflJlBer^in et al.| |2mj) must 
be required to c onvert atomic to molecular gas. Recently, 
iBlitz R.osolow skv (2004) note that the fraction of molec- 
ular gas may depend on a specific factor, such as the gas 
pressure, whilst the formation mechanism for GMCs is the 
same regardless of the molecular gas content. 

The concentration of GMCs to spiral arms can be ex- 
plained in one way by the increased local density ( Kennicutt 
Il989i ). Alternatively spiral density waves may take a 
more important role and lead to the direct trigge ring of 
star formati on (iRobertsI [l969l: iBonnell et al.1 l200d) Sev- 
eral models fBergin et al.' l2004l: iHollenbach McKeeiri979 ^: 
iKovama & Inutsuka 2000|) propose formation of molecu- 
lar hydrogen through shock compression of the ISM. Thus 
shocks from spiral density waves could account for the pref- 
erence of molecular gas in the spiral arms and the signifi- 
cance of spiral arms as sites of star formation. 

The properties of molecular clouds include significant 
structure on all observed scales an d supersonic bulk motions 
that follow t he relati on a., oc B^^'^ (|Mveralll983l : rDame et al.1 
ll98dlBrunJ 2003: El megreen Scalol2004l) . These kinemat- 
ics are generally taken to be turbulence and thought to gen- 
erate the internal structures. An alternative explanation is 
that structure in the interstellar medium generates the inter- 
nal velocity dispersion as it passes through a clumpy spiral 
sho ck feonnell et al. 2006). 

IWada et al.> ( 2002) have investigated self- gravitating 
galactic disks, suggesting that self- gravity of the disk gen- 
erates turbulence in the ISM. The effect of spiral shocks on 
the structure of non self- gravitating disks has also been de- 
scribed (Wada & Koda 2004; Dobbs & Bonnell 2006). Re- 
cent simulations of gravitational collapse in disk galax- 
ies have investigated the dependence of the star for- 
mation rate on the Toomre instability parameter and 
a Schmid t law dependence on the local surface density 
(iLi et al.1 |2005a). Molecular gas formation has generally 
been limited to smaller scale studies, e.g. cloud collapse 
(iMona pjhMLJ^Xajtanzid [1^911) and colliding or turbulent 
flows (jKovama Inutsukal200nHAudit Henneb elle 2005). 

In this paper, we look at the dynamics of a non-self- 
gravitating galactic disk subject to a spiral density wave. 
We show the formation of cl umpy molecular clou d struc- 
tures in the spiral arms. From lBergin et all (120041) we take 
a simple equation for the conversion of atomic to molecular 
gas which is applied to our galactic simulations. We identify 
the location of molecular hydrogen in the galactic disk and 
show the densities required for molecular gas formation. 



2 CALCULATIONS 

We use the 3D smoothed particl e hydrody namics (SPH) 
code based on the version by Benz ( Benzll99 0). The smooth- 



ing length is allowed to vary with space and time, with the 
constraint that the typical number of neighbours for each 
particle is kept near Nmigh ~ 50. Artificial viscosity is in- 
cluded with the standard parameters a = 1 and /3 = 2 
dMonaghan Lattanzioll 19851 : iMonaghanll 19921) . 



2.1 Flow through galactic potential 

The galactic potential includes components for the disk, 
dark matter halo and the spiral density pattern. The spi- 
ral potential is provided solely by the stellar background, 
which contains considerably more mass than present in the 
gas. Thus we ignore the gravitational feedback from the gas 
to the stellar disk. However, as the background potential is 
time dependent, the transfer of angular momentum between 
the potential and the gas is possible (Section 3.6). 
We represent the disk by a logarithmic potential 
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e.g. "Binnev & Tremaine^ ()l987l) that provides a fiat rotation 
curve with vo — 220 km s~^. The core halo radius is Rc = 
1 kpc and ^ = 0.7 is a measure of the disk scale height. A 
further component is the potential for the outer dark matter 
halo, solved from the density distribution 

Ph 



Phalo{r) 



1 + (r/r/,)2 



(2) 



where ph = 1.37x 10 ^ Mr:) pc ^ is the halo density and r/j = 
7.8 kpc the halo radius dCaldwell OstrikerlFaxl. For the 

spiral density pattern, we use a potential from lCox Gomez! 
(2002) 
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C7(3) = 8/157r 



(3) 

The number of arms is given by N and ro — 8 kpc, 
Rs = 7 kpc and H = 0.18 kpc are radial parameters. The 
pitch angle is a, the amplitude of the perturbation po and 
the pattern speed Qp. We take N = 4 for a 4 armed po- 
tential, and a pitch angle of a = 15^, comparable with the 
Milky Way (^Vallee 2005). The 4 armed potential also allows 
for more frequent passages of the gas through the spiral 
arms. The pattern speed, 2 x 10~^ rad yr~^ leads to a co- 
rotation radius of 11 kpc. The amplitude of the potential is 
^ 200 km^ s~^, with po = 1 atom cm"*^. This corresponds 
to a perturbation of 3 per cent to the disk potential, and an 
increase of 20 km s~^ as gas falls from the top to the bot- 
tom of the potential well. Overall the disk is in equilibrium, 
as the rotational velocities of gas in the disk balances the 
centrifugal force from the potential. The disk is in vertical 
equilibrium supported by an initial velocity dispersion. The 
resulting scale height of the disk is maintained throughout 
the simulation. 
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This paper only considers how hydro dynamic forces and 
galactic potential influence the flow. In particular we do not 
include the effect of self -gravity, magnetic fields and feed- 
back from star formation. The processes of heating and cool- 
ing of the gas are simplified by taking an isothermal equation 
of state. This approach is taken as our main objective is to 
assess the dynamics of spiral shocks rather than comprehen- 
sively model the ISM. 

2.2 Initial conditions 

We initially place gas particles within a region of radius 
5 kpc< r <10 kpc. The disk also has a scale height z ~ 
100 pc. We allocate positions and velocities determined from 
a 2D test particle run. The particles in the test particle run 
are initially distributed uniformly with circular velocities. 
They evolve for a couple of orbits subject to the galactic po- 
tential, to give a spiral density pattern with particles settled 
into their perturbed orbits. 

In the SPH initial conditions, we give the particles ve- 
locities in the z direction from a random Gaussian distribu- 
tion of 2.5% of the orbital speed. We add the same magni- 
tude velocity dispersion in the plane of the disk. The total 
mass of the disk is nominally 5 x 10^ M©, but as self-gravity 
is not included, the total gas mass can be scaled to higher 
or lower masses. The gas is distributed uniformly on large 
scales with an average surface density of S ^ 2 M© pc~^ 
and average density of 10~^ M© pc~^. For the higher res- 
olution run, most gas in the interarm regions has densities 
ranging from 10~^ M© pc~^ to 10~^ M© pc~^ whilst the 
initial peak density in a spiral arm is ~ 0.5 M© pc-^ (Fig- 
ure 1). The nominal surface density 2 M© pc~^) is lower 
than the gas surfac e density of 5 M© pc~^ suggested by 
dWolfire et al.l l2003'l for a radius of 8.5 kpc, although we 
mainly focus on modelling just the cold phase of the atomic 
gas. In Section 3.4 we consider the formation of molecular 
gas with different disk masses and therefore surface densi- 
ties. 

The number of particles is either 10^ or 4 x 10^. The 
higher resolution work was carried out on UKAFF (UK As- 
trophysical Fluids Facility) whilst the 10^ particle simula- 
tions were run on local computing facilities. We ran sim- 
ulations for between 200 and 300 Myr, where the orbital 
periods at 5 and 10 kpc are 150 Myr and 300 Myr respec- 
tively. The gas passes through multiple spiral shocks in this 
time and reaches a steady state. No boundary conditions 
are applied to our calculations. With 4 x 10^ particles, the 
disk is fairly well resolved in the spiral arms, where the ra- 
tio of the smoothing length (h) to the disk scale height (H) 
is approximately 0.4 (the disk scale height calculated as the 
height of the disk which contains 2/3 of the local disk mass). 
However the ratio h/H increases up to ~ 1 in the interarm 
regions. 

All calculations are isothermal with temperatures of 
10, 50, 10^, 10^, or 10^ K. Most previous simulations 
and analysis of spira l galaxies (e.g. I Wad a^ fcjCodal l2004l 
iKim OstrikeJ l2002l : |Dwarkadas & Balbus 1996) assume 
a sound speed of approximately 10 km s~^ corresponding 
to the hotter (10^ K) component of the ISM. However, for 
GMCs to form, gas entering the spiral arms must be cold 
ato mic clouds ([ Elmegre enI l2002l) or pre-existing molecular 
gas (jPringle et al.ll200lh . Since our aim of this work is to 
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Figure 1. The distribution of mass as a function of gas density at 
the beginning (diamonds), after 100 Myr (crosses) and at the end 
(280 Myr) (squares) of the 50 K simulation. The nominal density 
for molecular clouds is 2.5 M© pc~^. 

investigate molecular cloud formation in shocks, lower ISM 
temperatures are adopted for most of these calculations. We 
do not include self gravity, but estimate the Toomre param- 
eter for the disk to be Q ~ 2.7, and so the disk is sta- 
ble against gravitational collapse. In any case, the gas will 
reach a spiral arm before gravitational collapse can occur 
e.g. (BonneU et al. 2006). We refer to Li et al. (2005b) for a 
discussion of the stability of disks for different Q. 



3 DISC DYNAMICS AND CLOUD 
FORMATION 

We first describe briefly the overall structure for the differ- 
ent simulations performed. Section 3.1 then examines the 
formation of structure in the spiral arms, before we discuss 
the velocity dispersion across the disk and the distribution 
of molecular gas. 

Initially, the disk is approximately uniform with a 
smooth spiral perturbation. The widths of the spiral arms 
are even along the length of each arm, and the densities 
across the arms are fairly uniform. As the shock devel- 
ops, the density increases in the spiral arms. In the simu- 
lations with cold gas, substructure develops along the spi- 
ral arms (Figure 2). The gas in the arms becomes clumpy, 
and the arm width fluctuates. Densities of gas in the arms 
varies over approximately two orders of magnitude. The spi- 
ral arms become more ragged, and their width increases. 
Meanwhile dense clumps leaving the arms are sheared, lead- 
ing to the formation of spurs perpendicular to the arms 
(|Dobbs Bonnelll 12006'). The distribution of gas densities 
reaches an approximate equilibrium after 100 Myr (Fig- 
ure 1) with gas densities in the spiral arms reaching > 10 
M© pc~^. In contrast, the higher temperature simulations 
show broader, less dense spiral arms (see Dobbs & Bonneljj 
(2006)). The spiral shocks are much weaker for the higher 
temperature gas, since the change in velocity induced by the 
spiral potential is more comparable with the sound speed of 
the gas. For the 10^ K case, the sound speed is ^ 10 kms~^ 
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whilst the change in velocity across the spiral potential is 
^ 20 kms"^ (Section 2.1). The 1000 K simulation shows 
regular spaced more diffuse clumps along the spiral arms, 
but the 10^ K run is relatively featureless, with no obvious 
structure in the spiral arms. The difference in structure is 
due to the different strength of the shock and the spacing of 
the clumps i s determined by the spir al arm dynamics (Sec- 
tion 3.1 and lDobbs BonnelJ (l200(j) ). 
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Figure 2. Column density plots (g cm~^) showing a 5 kpc by 
5 kpc section of the disk (with Cartesian coordinates -9 kpc < 
X < -4 kpc and -5 kpc < y < kpc. The xy-coordinate grid is 
centred on the midpoint of the disk and remains fixed with time. 
Gas is flowing clockwise across the disk, i.e. from bottom to top. 
The time corresponding to each plot is i) 50 Myr, ii) 100 Myr 
and iii) 150 Myr. This simulation used 4 million particles with a 
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3.1 Structure formation in the spiral arms 

We associate the densest regions of our simulation (Figure 2) 
with molecular clouds, as they attain the required densities 
and, as we shall see in Section 4, have appropriate phys- 
ical conditions and lifetimes that they should be primarily 
molecular gas. The formation of dense molecular cloud struc- 
tures can be understood as being caused by the dynamics 
as interstellar gas passes through a spiral shock. Initially, 
the pre-shock gas contains small-scale structure due to the 
particle nature of SPH. The spiral shocks amplify any pre- 
existing structures such that the dumpiness of the gas in- 
creases with each passage through a spiral arm. This ampli- 
fication of dumpiness corresponds to changes in the angular 
momentum phase space density of the gas. The pre-shock 
phase space density is relatively smooth whilst the shocked 
gas displays much more structure (Figure 3). We provide an 
explanation of the generation of clumpy structure in terms 
of this change in the angular momentum of the gas each 
time it passes through a spiral shock. 

As the arms are regions of strongly enhanced density, 
it is evident that gas from a range of radii must interact 
and shock at any particular point in a spiral arm. These 
shocks disturb the initially smooth phase-space density of 
the gas particles by modifying the particles' angular mo- 
menta. In Figure 4 we illustrate this by plotting the evolu- 
tion of Lz, the 2;-component of the angular momentum, for 
a select number of particles as a function of time. Initially 
the particles are at the same azimuth, but are at a range of 
radii, with of course larger radii corresponding to larger Lz- 
We see that particles with higher initial values of L^, and 
therefore at larger radii, enter the shock earlier. This is be- 
cause the spiral arms are trailing and the particles are inside 
the co-rotation radius. As a particle enters such a shock it 
is typically at the outer radial range of its epicyclic motion 
wit hin the sp i ral po tential (see, for example. Figures 3 and 
9 in iRoberta ()l969h '). and so has lower angular momentum 
than the average for that radius. Thus when it enters the 
shock, it mixes immediately with material which has higher 
angular momentum. In Figure 4 the particles entry into the 
shock is marked by a sudden jump in the value of Lz. In 
Figure 5, we show that the magnitude of the jump is cor- 
related with the post-shock density the particle encounters. 
Thus the stronger the mixing, the bigger the jump. Each 
particle then travels along nearly paralle l to the shock for 
some distance, before it leaves the shock (|Roberti[l969|). As 
it moves along the shock, it stays roughly in phase with the 
rotating spiral potential, and is subject to a torque from that 
potential. Because these particles are inside co-rotation, the 
phasing of the arm relative to the potential ensures that the 
torque is a retarding one, and thus leads to a decrease in the 
value of Lz. Once the particle leaves one spiral arm it then 
proceeds at more or less constant Lz until it meets the next 
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Figure 3. The top figure shows particles in a section of the disk 
in coordinates (r, 0). is measured in radians clockwise round the 
disk, and thus the section is taken from the top of the disk. The 
bottom figure shows the corresponding z component of the angu- 
lar momentum plotted against 0. The evolution of the structure 
in the angular momentum phase space is evident in the bottom 
plot. The clumpy region shows the particles in the spiral shock. 
Gas below this has yet to pass through the shock and has a rel- 
atively smooth distribution. The post-shock gas lies above the 
clumpy region. The simulation has evolved for 50 Myr. 



one. As we would expect for particles inside co-rotation, we 
can see from Figure 4 that the decrease in Lz , which occurs 
as the particle slowly leaves the spiral arm, exceeds the in- 
crease it acquires on entry to the arm. Thus the net effect of 
the spiral potential, coupled with dissipation in the arms, is 
to produce a steady transfer of angular momentum from the 
particles to the potential. This steady trend of reducing the 
angular momentum of the particles is illustrated in Figure 6. 

We note that in the simulations with higher values for 
the temperature of the ISM, the shocks are weaker, and the 
distribution in Lz stays much more uniform. This is because, 
after a weaker shock the particles velocity is less affected and 
so the particle stays in phase with the rotating potential, 
and so subject to a retarding torque, for a smaller amount 
of time. 

Although the pre-shock ISM is inhomogeneous, we find 
that this dumpiness is amplified within the spiral arms. It 
is this amplification which gives rise to the structures within 
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Figure 4. The angular momentum (z-component) with time of 
10 selected particles with 6.8 kpc < r < 7.8 kpc. The particles 
are initially all at the same azimuth, but at different radii. 
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Figure 5. The jump is 2;-component of angular momentum 
(Lz(2) — Lz(l)) is plotted as a function of post-shock density 
(p(2)) for a number of particles. The pre-shock densities of the 
particles (not shown) all fall below the dotted line. There is a 
strong correlation between density enhancement in the shock and 
the magnitude of the jump in Lz- 



the arm (Figure 2), particularly in the molecular gas, as will 
be shown in Section 4.2. To see this in more detail. Figure 3 
(upper) shows the distribution of particles in a small sector 
of our simulation. Since is measured clockwise round the 
disk, particles are located in the top part of the disk, and are 
selected to cover a spiral arm. Approaching the spiral arm 
from below can be seen apparent waves of material. These 
are in fact the remnants of clumps formed in the previous 
arm, by now strongly elongated by the differential shear in 
the inter-arm region (See Figure 2). These inhomogeneities, 
which are transmitted from the previous arm, provide input 
for new inhomogeneities in the next arm. In Figure 3 (lower) 
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Figure 6. The number of particles (N(L2)) which have the value 
of the z component of the angular momentum . The total num- 
ber of particles for the whole disk is included at the beginning of 
the simulation (solid line), after 140 Myr (dashed line) and after 
280 Myr (dotted line). 



we plot the value of Lz for each of the particles in Figure 3 
(upper) as a function of the ^-coordinate. The region below 
the arm in the upper Figure 3 now spans a relatively small 
range in L^, and the phase-space density here is relatively 
smooth. Just above this is a region where the phase-space 
density drops. This corresponds to the rapid jump in L^, 
discussed above, which occurs as particles enter the spiral 
shock. Above this we see a patchy structure of phase-space 
density enhancements. It is clear that the inhomogeneity 
which occurs here, within the arm, is much greater than 
that seen in the pre-shock gas. 

The physical picture we have of what is happening here 
is as follows. We have seen that the gas, on entering the 
arm and being shocked, moves along the arm for a while 
(Roberts, 1969). It eventually leaves it to enter the inter-arm 
region, when its angular momentum is too high to follow the 
spiral arm to smaller radii. Imagine, then, as a simple model 
for the gas flowing along the arm, a set of masses flowing 
in a line along a smooth inclined tube under gravity. If the 
gas were all homogeneous, we would model this by adding 
masses to the tube at equal velocities, at a rate which is uni- 
form in time and space. To mimic material leaving the arm, 
we remove masses from the tube at positions which depend 
on their velocity (for example, we could remove particles 
when their velocities exceed some particular value Vmax)- In 
this case, at any one time the density of masses along the 
tube would remain smooth and independent of time. How- 
ever, to model an inhomogeneous pre-shock ISM, we would 
add masses to the tube in a non-uniform manner in space, 
and with a range of velocities. In this case, as the masses 
move along the tube they would tend (assuming the en- 
counters are dissipative) to bunch up, like traffic along a 
single lane highway. This bunching in velocity space can be 
seen in Figure 3 (lower). Since as before, we would remove 
masses from the tube at positions which depend on their 
velocities, it is clear that in this case, masses would not only 
tend to form clumps within the tube, but would also tend to 
be removed from the tube in clumps. As these clumps leave 
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Figure 7. The velocity dispersion plotted against azimuth (^, 
radians) for a ring whereby r=7.5 kpc and Ar=200 pc. The figure 
shows the dispersion of the v_\_ = -^J {v^ ~^ ^y) (solid) and the 
dispersion of Vz (dotted). (0 is measured clockwise round the disk, 
i.e. in the direction of gas fiow). Insufficient resolution is likely to 
produce damping in the interarm regions. The figure corresponds 
to a time of 100 Myr in the simulation. 

the spiral arms, they are sheared out to form the spurs and 
feathering seen in Figure 2. 

Using these ideas we can obtain a rough estimate of 
the linear scale of the inhomogeneities we expect to find in 
the spiral arms. The typical time a particle spends within a 
spiral arm can be estimated from Figure 4 as being approx- 
imately tarm ~ 6 X 10^ ycars. The typical velocity spread 
of particles within a spiral arm can be seen from Figure 7 
to be approximately v± ^ 10 km s~^. Then the average 
scale-length of density structure generated by the piling up 
process would be L tarmV± ^ 600 pc. This agrees roughly 
with the mean distance between the high-density peaks, 
most easily viewed from the molecular gas density (yellow) 
along the arm in Figure 9(iv). 

3.2 Velocity dispersion 

In addition to introducing structure in the phase space den- 
sity of the gas, spiral sho cks can also induce large veloc- 
ity dispersions in the gas ('Ron nell et Zll2nn6^ . We deter- 
mine the 2D velocity dispersion in the plane of the disk as 
a function of azimuthal angle. This is calculated by taking 
a 200 pc width ring centred at r = 7.5 kpc. We divide this 
ring into 100 segments, each of length rAO = 50 pc and area 

9 X 10^ pc^ and calculate the velocity dispersion of the gas 
in each segment. Figure 7 shows the variation of the velocity 
dispersion for the 50 K simulation. The velocity dispersion is 
calculated after 100 Myr, for the velocity component in the 
plane of the galaxy, v± = ^ (v% + t'y), and the z component 
of the velocity. There is a large increase of 10 km s~^ in the 
dispersion of v±_ where the spiral arms occur. In comparison 
a velocity dispersion of 7-1 km s~^ is observed for molecu- 
lar clouds of sizes 20- 5 pc JSolomon et al.ll98^lDame et al.l 
1986). Bo nnell et al.l (|2Q06h find a velocitv dispersion up to 

10 km s~^ in simulations of gas passing through a spiral 
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shock. This velocity dispersion is due to random chaotic in- 
ternal motions induced in clumpy gas subject to a shock. 
The velocity dispersion quickly decreases upon leaving the 
arms. A further decrease in the velocity dispersion occurs 
in the interarm regions. The ratio in the smoothing length 
compared to the disk scale height (h/H) approaches 1 in 
the interarm regions. This potentially leads to damping of 
velocities due to insufficient resolution, particularly in the 
Vz velocity component. 



4 MOLECULAR GAS 

We have seen above that spiral shocks produce highly struc- 
tured and dense gas that can be associated with molecular 
clouds in spiral arms. In addition, we would like to know 
whether such a formation mechanism is able to account for 
the formation of molecular gas from initially atomic gas or 
if such gas would need to be primarily molecular before en- 
tering the shock. In order to do this, we perform a rough 
estimate of the formation (and destruction) of molecular 
gas as the gas enters and leaves the spiral arms. We post- 
process the simulation by using the local gas properties de- 
termined by the isothermal simulations and calculate the 
time dependent formation rate of H2 on dust grains as well 
as the dissociation rates by UV radiation and cosmic rays. 
At each time step from the simulation, the formation and 
destruction rate is calculated for each gas parcel (particle) 
and added to its previous molecular fraction. In this way, 
we can estimate the spatial and temporal evolution of the 
molecular gas throughout the simulation. 

Since the gas in our models is isothermal, the forma- 
tion of molecular hydrogen is not treated self-consistently. 
In particular there is no heating o r cooling;, as de- 
scribed in previous smaller scale analvsis (|B erp;in et al.l2004l : 
lAudit Hennebelle"2nn5': 'Kovama Inutsuka"200n\ As a 
first attempt to analyze the overall distribution of molecular 
hydrogen, and the densities required for H2 formation, we 
post-process our results to find the time-dependent molecu- 
lar density of gas in the disk. This is a simplistic approach 
but identifies the most dense gas in our simulations where 
molecular hydrogen forms. The reader should bear in mind 
however that heating from feedback and supernovae explo- 
sions will have some effect on these results. 

We use a simple algorithm for the evolution of molecu- 
lar hydrogen (|Bergin et al.ll2004l) determined for molecular 
cloud formation behind shock waves. The rate of change of 
H2 density is given by 



dn{H2 
dt 



= Rgr{T)npn{H) - [Car + Cd^ss{N{H2),Av)HH2) 

(4) 

where n is the number density (of atomic or molecular hy- 
drogen), N is the column density (of atomic or molecular 
hydrogen) and T the temperature. The total number den- 
sity is rip, where rip = n{H) + 2n(H2). Molecular hydrogen 
is produced through condensation onto grains, determined 
by the first term on the right of eqn(4). Rgr is the forma- 
tion rate on grains, which is proportional to T°"^. The effi- 
ciency for the formation rate on grains is 5* = 0.3 , although 
recent results suggest the efficien cy is hig!;her than this fo r 
temperatures of less than 90 K (jCazaux TielenslEool . 
Destruction of H2 occurs through cosmic ray ionization and 



photo-dissociation, represented by the second term on the 
right of eqn(4). The cosmic ionization rate, Ccr is assumed 

constant (Ccr = 6 x 10"^^ s' M- 

Follo wing iBergin et al.l feOOA. we adopt the approxi- 
mation of iDraine BertoldJ (Il99^ for the H2 dissociation 



rate: 



Cd^ss{N{H2)) = Cd^ssfsh^eld{N{H2))e 



'Cd^ss{0) 



fsMeld{N{H2)) 



0.965 



+ ■ 



0.035 



(l + x/3)2 (l + x)0-5 
X exp[-8.5 X 10"^(l + x)°-^] 



(5) 



where x = N{H2)/5 x lO^^cm"^ 



The function f shield approximates the rate of photodissocia- 
tion summed over all energy levels. Self shielding is included 
through the dependence on column density and Td,iooo is 
the dust extinction at 1000 A. The constant Cd^ss(0) = 
4.17 X 10~^^ s~^ is a measure of the strength of the UV 
field. This formulation assumes that the background UV 
spectrum is the same as that of a BO star, and that the 
radiation flux is constant everywhere. The strength of the 
flux is measured by the magnitude of Cdiss{0). We note that 
this is a very approximate measure of the the photodissoci- 
ation rate, which should ideally vary depending on the local 
level of star formation. Furthermore, we only produce an es- 
timate of the column density N(H2) in our calculations by 
multiplying the local density by the scale height of the disk. 

We divide Eqn(4) by Up and calculate the change in the 
molecular gas fraction, d{n{H2)/np), for each particle after 
each timestep (2 Myr). Assuming n{H2) = everywhere 
initially, we add the change in the molecular gas fraction to 
the ratio n(H2)/np determined for the previous timestep. 
At each time frame, this ratio in the number densities can 
be converted to a volume density of molecular hydrogen. As 
shown in Figure 8, Eqn(4) mainly selects the most dense 
particles, although the molecular gas fraction is also time 
dependent. Since the formation of H2 is not instantaneous, 
the molecular gas fraction increases with the time gas spends 
in the denser spiral arms. 



4.1 Molecular gas formation at different 
temperatures 

We find that H2 formation in our simulations requires that 
the (isothermal) gas be cold, with temperatures of order 
100 K or less. Of order 10% of the gas becomes molecular 
with these temperatures whereas when the gas is 1000 K 
the molecular fraction drops to around 1% and is virtually 
zero at 10"^ K. We can thus conclude that the spiral shock 
formation of H2 is likely when the ISM is cold. The fraction 
of molecular gas was found to decrease with lower resolu- 
tion. For the 50 K highest resolution simulation, the total 
molecular gas in the disk reaches a maximum of ^ 13%. The 
fraction of molecular gas was found to decrease with lower 
resolution. For 1 million particles, the amount of molecular 
hydrogen formed was very similar in the 10, 50 and 100 K 
cases, reaching a maximum of ~ 8% of the total mass. 

The fact that the gas is assumed to be isothermal pre- 
cludes a detailed study of the temperature dependency on 
H2 formation as in reality this would depend on the details 
of the equation of state. For example, the efficiency of grain 
formation, assumed here to be a constant S=0.3, can be as 
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Figure 8. The distribution of the molecular fraction of particles 
is shown versus their density for the 50 K simulation, after 100 
Myr. 1/lOth of particles are plotted, i.e. 4 x 10^. There is a clear 
cut off at Up ~ 10 g cm""^, below which the gas is totally 
atomic. The large spread in the H2 fraction for a given density 
represents the time evolution of the H2 fraction from when low 
density gas first enters the spiral shock. 

high as S=l for T < 20 a nd as low as 5' - 0.1 for T > 100 K 
(|Cazaux Tielenl 12004) . What we can tell is that hot gas 
(T ^ 1000 K) does not form significant H2 due to the shock 
dynamics. 

When the gas is hot, the shock produced by the spiral 
potential is weaker, so the density of gas in the spiral arms 
is reduced. The density in the lO'^ K simulation does not 
increase above 10~^^ g cm""^ which is too low for molecular 
hydrogen to form (Figure 8). The density of gas in the spiral 
arms is insufficient for self-shielding to occur and molecu- 
lar gas is r eadily phot odissociated. Similar conclusions were 
reached bv lBergin et a l. ( 2004), who require pressures a few 
times higher than the average interstellar values to allow ef- 
fective self shielding. We do not inclu de cooling, although 
for st rong shocks performed at 10^ K ((Kovama Inutsukal 
12000!) only a thin layer of molecular hydrogen forms. There- 
fore at higher temperatures, it is unlikely that spiral shocks 
could produce densities sufficient for significant molecular 
hydrogen formation. 

4.2 High resolution simulation 

In the remainder of the paper, we focus on the 50 K highest 
resolution simulation. We show the location of H2 when the 
ISM temperature is 50 K in Figure 9. There are 4 x 10^ 
particles giving a mass resolution of 125 per particle. 
This figure is taken after 100 Myr, and shows the global 
distribution of molecular gas, and the small scale structure 
of molecular clouds. 

Molecular hydrogen starts to form immediately in the 
spiral arms, the most dense parts of the disk. H2 is first 



situated in a narrow strip, where the shock from the spiral 
density wave occurs. The density of gas leaving the arms is 
insufficient to prevent molecular gas quickly becoming pho- 
todissociated, hence the molecular gas is solely confined to 
the shock in the spiral arm. This strip of molecular hydrogen 
is smooth and continuous along the spiral arm. However, as 
the density of the shock increases, denser structures become 
apparent in the arms. These break up into separate clumps 
of molecular gas, which we identify as 'molecular clouds'. 
These molecular clouds are produced simply by the kine- 
matics of the gas flowing in the spiral potential as described 
in Section 3.1, as self gravity is not included. As the sim- 
ulation progresses, dense structures leaving the spiral arms 
become sheared. Some regions of molecular gas leaving the 
spiral arms are then dense enough for self shielding to re- 
duce the rate of photodissociation. These molecular clouds 
survive further into the inter arm regions. 

4.3 Molecular gas and spiral structure 

Figure 9 shows that molecular gas is largely conflned to the 
denser spiral arms. We plot the average molecular and total 
(molecular and atomic) gas densities with azimuth in Fig- 
ure 10, at 2 different times. As in Section 3.2, the densities 
are calculated from a ring centred at r = 7.5 kpc divided 
into 100 segments. This plot only displays the mean molec- 
ular and atomic densities - the peak densities are 10-20 times 
higher. The interarm densities of H2 are 10^ times smaller 
than in the spiral arms after 100 Myr, whilst the total den- 
sities are ^ 50 times smaller. After 100 Myr, 10% of the 
gas in the ring is molecular. Approximately 70% of the total 
amount of gas is situated in the spiral arms, of which 20% is 
molecular. Figure 10 suggests that the existence of interarm 
molecular gas is unlikely, although we examine the effect of 
changing the disk mass and photodissociation rate next. 



4-. 3.1 Dependence of H2 formation on disk mass and 
photodissociation rate 

The mass of H2 in eqn (1) is determined by the density, pho- 
todissociation rate and temperature. The cosmic ray ioniza- 
tion rate, ^cr = 6 x 10~^^ s~^, is found to have a minimal 
effect on the fraction of molecular to atomic gas. Since the 
H2 mass is determined post process, calculations can be re- 
peated with a different disk mass (as there is no self grav- 
ity) and/or different photodissociation rate. As expected, 
increasing the mass or decreasing the photodissociation rate 
considerably increase the formation of H2. The mass of H2 
increases considerably to 1.4 x 10^ M©, 28% of the total gas, 
when the disk mass increases by a factor of 10. In compar- 
ing with the photodissociation rate, the same approximation 
(eqn (4)) was applied, but increased or decreased by a con- 
stant factor. When the photodissociation rate decreases by a 
factor of 10, the molecular mass increases to 1.25 x 10^ M0, 
25% of the total gas mass. 

The disk mass has a substantial effect on the arm- 
interarm fractional abundance of H2. When the mass of the 
disk increases to 5 x 10^ M©, nearly all the gas in the spiral 
arms is molecular at the end of the simulation (Figure 11). 
Whilst the interarm ratio of the total density to H2 den- 
sity also increases, only about 1/lOth of the interarm gas is 
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Figure 9. Column density plots (g cm~^) showing density of molecular hydrogen (red) against overall density (black and white). Length- 
scales of plots are i) 20 kpc x 20 kpc, ii) 10 kpc x 10 kpc, iii) 5 kpc x 5 kpc (-8.5 kpc < x < -3.5 kpc, -8.5 kpc < y < -3.5 kpc), iv) 3 
kpc X 3 kpc (-6.5 kpc < x < -3.5 kpc, -6.5 kpc < y < -3.5 kpc). The xy-coordinates assume a Cartesian grid centred on the midpoint of 
the disk which remains fixed with time. Gas is flowing clockwise across the disk. 



molecular. The increase in both the arm and interarm re- 
gions is due to higher densities increasing the formation of 

H2. 



is arriving at the spiral arm which has not fully dissociated 
from the previous spiral arm passage. 



The photodissociation rate has a more pronounced ef- 
fect on the interarm ratios. When the photodissociation rate 
is a tenth of that taken from lPraine Bertoldil (jlQQGl) . ap- 
proximately a quarter of the interarm gas is molecular at the 
end of the simulation (Figure 11). By contrast the amount 
of H2 in the spiral arms changes much less. This difference 
at a lower photodissociation rate arises when the gas is not 
fully dissociated between the spiral arms. Molecular gas can 
survive for longer into the interarm regions leading to some 
molecular structures between the arms. This is illustrated in 
Figure 12, which displays part of the 50 K simulation where 
the molecular gas density is calculated with the lower pho- 
todissociation rate. Clumps of molecular gas are shearing 
away from the spiral arm. In addition, molecular material 



A more complete treatment of H2 formation will give a 
more accu rate idea of the quantitative analysis of molecular 
hydrogen (|Glover Mac Lowl l2QQ6al lb|) . Including heating 
may reduce the amount of molecular gas, but we defer com- 
menting on the effects of feedback, and varying the photodis- 
sociation rate to depend on star formation to future studies. 
These results indicate nevertheless that the molecular gas 
is strongly confined to the spiral arms. The possibility that 
molecular hydrogen is transmitted between spiral arms can- 
not be ruled out though, e.g. if the photodissociation rate 
is sufficiently low. The effects of increasing the disk mass or 
decreasing the level of photodissociation are likely to hold 
when including heating and cooling of the gas. 
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Figure 10. Density of molecular gas (solid) and total density 
(dotted) are plotted against azimuth (^), where (in radians) is 
the angle measured in the direction of flow. The density plotted is 
the average over an annulus centred at r=7.5 kpc with 5r=200 pc. 
The time of the simulation is 100 Myr (top) and 290 Myr (bot- 
tom). 



5 OBSERVATIONAL COMPARISON OF 
MOLECULAR CLOUD PROPERTIES 

One of the objectives of these simulations is to see how well 
the spiral shocks reproduce the observations of molecular 
clouds. We can directly compare our results with molecu- 
lar clouds in the Milky Way. We consider the distribution 
of clouds from a local standard of rest comparable to the 
Sun, and a further LSR at the edge of the disk. Several 
surveys of molec ular clouds in the Milky Way have bee n 
undertaken (e.g. iDame et al.lll986l . '2001: W et aLlEonl . 
We compare our data with 2 s ets of results: i) th e FCRAO 
CO survey of the outer galaxy (jHever et al.ll200ll) . This sur- 
vey detected molecular clouds between the Galactic longi- 
tudes 102^^.49 and 141^.54 which are located predominantly 
in t he local and Perseus arms; ii) a CO molecular line sur- 
vey (lYang et al.lEo0 2'l. which detected cold infrared sources 
over the northern plane of the Galaxy (longitudes between 
0^ and 260^). 

Our calculations give the fraction of molecular gas for 
each particle in the galactic disk. We use a simple algorithm 
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Figure 11. Fractional abundance of H2 in the arm (solid) and 
interarm regions (dashed) compared to the total disk mass (top) 
and photodissoclation rate (bottom). These figures are taken from 
the end of the simulation. 



to determine where molecular clouds are located. We select a 
section of 1/8 of the disk as a sample, and divide this region 
into a 5 pc resolution grid. We calculate the total fraction 
of molecular gas in each cell and retain cells which contain 
over 10 M© of H2. We classify all groups of adjacent cells 
satisfying this threshold as a molecular cloud. The properties 
of each 'cloud', including mass, position and velocity are 
tabulated. We determine the Ylsr of each molecular cloud 
-the locations of each LSR we have chosen are shown on 
Figure 13 compared to the overall distribution of molecular 
clouds. 

In Figures 14 & 15, typical results from our calcula- 
tions are compared with the FCRAO CO survey and the CO 
molecular line survey (where the LSR is situated at Carte- 
sian coordinates (0,6.5)kpc). Our grid based method gave 
a total of 5044 objects within the sample criterion. Below 
these results corresponding data from the FCRAO CO sur- 
vey are shown, where there are 7724 objects and the molecu- 
lar line survey, where there are 1331 objects. The data from 
our model and the 2 surveys show clear bands, where the 
molecular clouds have preferential radial velocities due to 
the spiral arm structure. The molecular clouds in our model 
are more confined to the spiral arms than in the observa- 
tional data, particular for the outer arm. We have assumed 
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Figure 12. A 3 kpc by 3 kpc region of the disk showing the 
transfer of molecular hydrogen b etween spiral arms . The pho- 
todissociation rate is 1/10 that of iDraine BertoldJ and 
as a consequence, there is more molecular hydrogen surviving 
between the spiral arms. The column density of molecular hy- 
drogen (red) and total density (black and white) are displayed 
(same scale as Figure 4). The times of each plot are i)100 Myr, 
ii)120 Myr and iii)140 Myr. Plot i) is the same time, region and 
scale as Figure 4(iv), but using a lower photodissociation rate. 
The increased density of H2 along the arms leads to molecular 
structures within the interarm regions which are absent from the 
simulation in Figure 4. 




Figure 13. The position of each local standard of rest (crosses) 
is compared to the overall distribution of molecular clouds. The 
LSR situated at (0,6.5)kpc is comparable with the position of the 
Sun, in an interarm region in the mid-part of the disk. The second 
LSR (0,10.5)kpc produces a distribution of clouds as viewed from 
the edge of the disk. 

a corotation radius of 11 kpc, which may reduce the spiral 
arm shock at larger radii. 

Figure 16 displays the molecular clouds as viewed from 
outside the disk, at coordinates (0,10.5)kpc. This is more 
comparable with Stark & Lee (2005), who observe clouds 
in the inner regions of the Milky Way. However their sur- 
vey incorporates 5 spiral arms, whereas only 2 are clear in 
Figure 16. 



6 CONCLUSIONS 

The simulations in this paper demonstrate the formation of 
molecular clouds in spiral galaxies, and specifically the role 
of spiral shocks. The compression as the gas passes through 
a spiral shock produces high density regions that are appro- 
priate for the formation of molecular hydrogen. A necessary 
condition for this to occur is that the gas is cold (T < 100 K). 
Warm gas (T ^ 1000 K) produces little molecular hydrogen. 

The spiral shock also produces significant structure in 
the gas along the arms. Cold gas is unstable to clumping in 
the shocks due to the semi-random encounters with lower 
and higher angular momentum gas. These denser structures 
are predominantly molecular and therefore identifiable as 
molecular clouds. An increased velocity dispersion is appar- 
ent in the spiral arms and thus molecular clouds, of similar 
magnitude to observations. We associate this increase in ve- 
locity dispersion with the generation of chao tic motions in 
the ISM as gas passes through clumpy shocks teonn ell et al.l 
2006). We find only a small increase though in velocity dis- 
persion in the vertical direction. Our simulations are not 
well resolved in the z direction, and we would probably not 
expect to find hig h velocities in this direction witho ut stellar 
energy injection (Ide Avillez BreitschwerdtJl2005l) . 
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Figure 14. V^^^ plotted against galactic longitude for molecular 
clouds from these simulations and the FRCAO CO survey. The 
top plot shows the simulated data, with a LSR at coordinates 
(0,6.5)kpc. T he lower plot uses data from the FCRAO CO survey 
jHever et al.lll998l 



The total percentage of molecular hydrogen is typi- 
cally ~ 10% for gas less than 100 K and a total disk gas 
mass of 5 X 10^ M©, (corresponding to a surface density of 
~ 2 M0 pc~^), but this figure roughly doubles with a sur- 
face density of 20 Mq pc~^ or lower photodissociation rate. 
We note that these figures are potentially lower limits on 
the amount of molecular hydrogen formed in these models, 
since we find that the total percentage of molecular hydro- 
gen increases with resolution. A caveat is that these results 
do not include realistic heating and cooling, which are likely 
to have some effect. 

Since molecular gas formation is determined by den- 
sity, and the time spent in dense regions, molecular gas is 
largely confined to the spiral arms. Typically we find the ra- 
tio of arm to interarm molecular hydrogen is 2 or 3 orders of 
magnitude. However a reduced photodissociation rate allows 
the survival of molecular gas into the interarm regions, and 
a ratio of arm to interarm molecular hydrogen less than 10. 
This leads to the possibility of molecular gas passing from 
one arm to the other and potentially determining subsequent 
molecular cloud formation. 

The distribution of molecular clouds in our simulations 
is consistent with observational surveys of molecular clouds 
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Figure 15. V^^^ plotted against galactic longitude for molecular 
clouds from these simulations and the CO molecular line survey. 
The top plot shows the simulated data, with a LSR at coordinates 
(0,6.5)kpc. The lower plot uses data from I Yang et al.l (|200A . 
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Figure 16. Vrms plotted against galactic longitude from simu- 
lated data. The LSR was chosen at coordinates (0,10.5)kpc 
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in the Milky Way. The confinement of molecular clouds to 
the spiral arms is more comparable to the outer Milky Way, 
where the ratio of arm to interarm clouds is high. The resolu- 
tion is insufficient to gain accurate information about these 
clouds. Higher resolution simulations will determine whether 
the properties of individual molecular clouds (e.g. mass or 
internal velocity dispersion) can also be reproduced. 
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